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A  two-dimensional,  axisymmetric  transient  computational  fluid  dynamics  model  is  developed  for  an 
intermediate  temperature  micro-tubular  solid  oxide  fuel  cell  (SOFC),  which  incorporates  mass,  species, 
momentum,  energy,  ionic  and  electronic  charge  conservation.  In  our  model  we  also  take  into  account 
internal  current  leak  which  is  a  common  problem  with  ceria  based  electrolytes.  The  current  density 
response  of  the  SOFC  as  a  result  of  step  changes  in  voltage  is  investigated.  Time  scales  associated  with 
mass  transfer  and  heat  transfer  are  distinguished  in  our  analysis  while  discussing  the  effect  of  each 
phenomenon  on  the  overall  dynamic  response.  It  is  found  that  the  dynamic  response  is  controlled  by  the 
heat  transfer.  Dynamic  behavior  of  the  SOFC  as  a  result  of  failure  in  fuel  supply  is  also  investigated,  and 
it  is  found  that  the  external  current  drops  to  zero  in  less  than  1  s. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  are  among  the  prominent  candidates  for 
an  alternative  power  generation  method  due  to  fuel  flexibilities 
[1,2]  and  their  capabilities  to  be  used  in  cogeneration  applica¬ 
tions  for  increased  system  efficiencies  [3-5].  On  the  other  hand 
SOFCs  with  micro-tubular  design  are  accepted  as  more  favorable 
than  planar  SOFCs  because  of  their  better  sealing  capabilities  and 
enhanced  thermo-cycling  characteristics  [6].  Improvement  in  ther¬ 
mal  management  and  reduction  of  the  losses  [7,8]  due  to  their 
sub-millimeter  dimensions  are  also  reported  as  a  superiority  over 
the  planar  design. 

Although  significant  amount  of  work  exist  on  micro-tubular 
SOFCs  [9-16],  these  devices  are  not  employed  in  any  practical 
applications  yet  due  to  lack  of  operation  optimization  to  maintain 
electrochemical  performance  and  mechanical  durability  simulta¬ 
neously.  This  is  also  related  to  the  difficulty  in  experimentation  as 
the  data  such  as  temperature  and  current  distribution,  is  not  so 
easily  measurable  due  to  the  smaller  dimensions  and  elevated  oper¬ 
ating  temperatures.  Consequently,  modeling  becomes  the  primary 
tool  to  develop  essential  insight  into  operation  since  one  can  acquire 
information  at  a  specific  point  inside  the  cell  through  modeling. 
For  example,  it  is  valuable  to  predict  where  the  maximum  temper¬ 
ature  occurs  during  the  fuel  cell  operation  from  a  thermal  stress 
perspective.  Furthermore,  employing  modeling  to  understand  the 
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dynamic  behavior  of  the  cell  is  necessary  because  it  is  vital  to  main¬ 
tain  mechanical  stability  of  the  cell  during  any  proposed  operation 
hence  it  is  essential  to  track  the  changes  in  temperature  with  time 
at  any  location  inside  the  cell.  Helping  to  understand  the  tran¬ 
sient  physical  phenomena  inside  the  fuel  cell,  modeling  becomes 
an  important  tool  to  predict  the  response  of  fuel  cell  parameters  to 
any  load  change  and  then  can  be  used  to  suggest  a  control  scheme 
for  the  fuel  cell  operation.  Modeling  can  also  be  utilized  to  identify 
the  fuel  cell  behavior  at  some  critical  conditions  such  as  failures  in 
system  components. 

While  there  are  a  number  of  significant  modeling  efforts, 
both  in  steady  state  [17-24]  and  in  transient  [25-33],  transient 
SOFC  models  typically  are  not  as  elaborate  as  steady-state  mod¬ 
els.  Some  of  the  transient  models  are  lumped  models  [34,35] 
neglecting  ah  the  spatial  variations.  They  are  mainly  developed 
for  control  or  to  simulate  the  fuel  cell  as  part  of  the  larger  sys¬ 
tem.  Transient  models  incorporating  transport  phenomena  vary 
from  one-dimensional  to  three-dimensional  models.  3D  models 
are  developed  generally  for  planar  SOFC,  whereas  exploiting  the 
axial  symmetry  of  the  tubes,  two-dimensional  models  are  preferred 
and  most  of  the  time  sufficient  to  represent  a  SOFC  with  a  tubular 
design. 

Achenbach  presented  one  of  the  first  transient  studies  on  SOFCs 
[26].  He  developed  a  3D  transient  model  to  investigate  the  volt¬ 
age  responses  of  a  planar  SOFC  to  certain  load  changes.  Ioara  et  al. 
developed  a  1-D  model  for  a  planar  direct  internal  reforming  SOFC 
in  which  they  concluded  that  neglecting  the  variation  of  the  gas 
stream  properties  may  lead  to  incorrect  dynamic  response  predic¬ 
tions  [27]. 
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Nomenclature 

Symbols 

k  thermal  conductivity 

T  temperature 

Cp  specific  heat 

u  velocity  vector 

Q.  volumetric  heat  source 

p  pressure 

1  unit  vector 

g  gravity 

R  volumetric  consumption 

w  species  mass  fraction 

x  species  mole  fraction 

n  number  of  species 

F  Faraday  constant 

R  universal  gas  constant 

i  exchange  current  density 

A  pre-exponential  constant 

EA  Electrode  activation  energy 

kb  Boltzmann  constant 

h  enthalpy 

M  molecular  weight 

rp  average  pore  radius 

Greek 

p  density 

/x  dynamic  viscosity 

k  permeability 

o  conductivity 

cp  potential 

y  concentration  dependency 

a  transfer  coefficient 

r]  overpotential 

s  porosity 

r  tortuosity 

Subscripts 
a  anode 

c  cathode 

elc  electrolyte 

j  species 

k  species 

i  ionic 

e  electronic 

eq  equilibrium 

eff  effective 

Kn  Knudsen 

Superscripts 
for  formation 

in  inlet 

reac  reacted 

T  transpose 

Th  thermal 


Jia  et  al.  provides  detailed  analysis  of  the  effects  of  operating 
parameters  on  the  steady  state  and  transient  performance  of  a 
tubular  SOFC  [29].  Employing  a  1 D  model  they  represented  the  con¬ 
servation  laws  with  a  control  volume  approach.  In  a  recent  study  of 
Barzi  et  al.  [28]  dynamic  responses  of  a  tubular  SOFC  are  predicted 
during  the  startup.  Their  2D  model  incorporates  mass,  momentum 
and  species  balances  accompanied  with  a  circuit  representation  of 
charge  balances.  Bhattacharyya  et  al.  compared  the  dynamic  behav¬ 


ior  of  the  cell  with  the  experiments  in  their  2D  model  of  a  tubular 
SOFC  [30].  Along  with  investigating  dynamic  response  of  the  cell  to 
the  changes  in  voltage,  they  also  predicted  the  response  of  the  cell 
to  the  changes  in  hydrogen  flow  rate. 

Ota  et  al.  compared  transient  characteristics  of  a  standard  tubu¬ 
lar  cell  with  a  micro-tubular  cell  with  a  modeling  framework 
presented  therein  [31].  Their  modeling  framework  is  based  on 
simplifications  instead  of  taking  into  account  full  coupling  of  the 
sophisticated  transport  phenomena.  They  reported  that  time  scales 
of  a  standard  tubular  cell  to  a  specific  voltage  response  is  six  times 
larger  than  that  of  a  micro-tubular  cell.  Another  modeling  study 
on  micro-tubular  SOFCs  is  presented  by  Nehter  [32].  In  this  study, 
he  compared  a  common  micro-tubular  cell  with  a  cascaded  one. 
Although  localized  temperature  and  species  concentrations  are 
provided  in  his  2D  axial  symmetric  model,  momentum  balance  and 
multicomponent  species  transfer  are  not  included.  Mass  balance  is 
carried  out  in  a  simple  way  via  algebraic  equations  describing  the 
electrochemical  and  shift  reactions. 

Although  there  are  many  simplifications  in  the  models  of  Ota 
et  al.  and  Nehter,  their  studies  are  significant  since,  to  the  knowl¬ 
edge  of  the  authors,  they  are  the  only  modeling  efforts  emphasized 
dynamics  of  micro-tubular  SOFCs.  Since  the  characteristics  of  a 
micro-tubular  SOFC  are  very  different  than  a  standard  tubular  cell, 
there  is  still  need  for  a  rigorous  model  to  study  the  dynamics  of 
a  micro-tubular  cell.  Therefore  in  our  study,  we  extend  our  previ¬ 
ous  modeling  framework  [36,37]  to  a  transient  model  to  study  the 
dynamic  behavior  of  a  micro-tubular  SOFC.  The  model  incorporates 
mass,  species,  momentum  and  energy  balances  along  with  electro¬ 
chemical  kinetics  and  ionic  and  electronic  charge  conservations. 
Taking  into  account  full  coupling  between  these  physical  phenom¬ 
ena,  model  equations  are  solved  in  a  distributed  space  and  time 
domain.  Due  to  the  axial  symmetry  of  the  tubular  geometry,  the 
model  is  constructed  in  a  two-dimensional  axisymmetric  domain 
by  assuming  the  anode  and  the  cathode  current  collectors  are  uni¬ 
formly  distributed  on  the  electrode  surfaces.  Simulations  are  then 
carried  out  to  analyze  the  dynamic  response  of  the  micro-tubular 
fuel  cell  and  predict  the  evolution  of  the  temperature  and  hydrogen 
concentration  profiles  as  a  result  of  voltage  changes.  In  addition,  we 
investigate  a  failure  scenario  in  which  the  flow  in  the  fuel  channel 
is  cut  off  suddenly  due  to  a  probable  malfunction  in  the  fuel  supply 
system. 


2.  Mathematical  model 

As  shown  in  Fig.  1,  model  geometry  includes  the  anode  tube 
and  the  fuel  channel,  the  electrolyte,  the  cathode  and  the  air  cham¬ 
ber  and  the  supporting  alumina  tubes.  Fig.  1  (not  drawn  to  scale) 
shows  the  cross-section  of  the  axial  symmetric  geometry.  The  actual 
model  geometry  can  be  visualized  by  revolving  this  cross-section 
around  the  symmetry  axis,  i.e.  centerline  of  the  anode  tube.  Mass, 
species  and  momentum  conservation  equations  are  solved  in  the 
gas  channels  and  the  porous  electrodes,  and  energy  conservation 
is  applied  in  the  entire  domain.  Ionic  charge  balance  is  applied  in 
the  anode,  the  electrolyte  and  the  cathode  whereas  the  electronic 
charge  balance  is  applied  in  the  anode  and  cathode. 

In  this  study,  we  model  the  micro-tubular  cells  which 
are  fabricated  and  characterized  by  the  New  Energy  Devel¬ 
opment  Organization  (NEDO)  of  Japan  under  the  Advanced 
Ceramic  Reactor  Program.  The  anode  supported  cells  employs  a 
Ceo.gGdo.i  0!  95_5  (GDC)  electrolyte  coated  on  a  NiO-GDC  anode,  and 
Lao.8Sro.2Coo.6Feo.4O3  (LSCF)-GDC  cathode.  It  is  well-known  that 
GDC  electrolytes  are  prone  to  internal  current  leaks.  GDC  becomes 
reduced  and  electronically  conductive  under  certain  conditions, 
especially  at  elevated  temperatures  [38].  Fraction  of  the  electrons 
generated  at  the  anode  is  then  transferred  to  the  cathode  through 
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Equations 

1.  Continuity 

2.  Momentum 

3.  Species 

4.  Heat 

5.  Electronic 
Charge 

6.  Ionic  Charge 
*  Momentum 

equation  in 
porous  media 


Fig.  1.  Model  geometry  and  the  equations  solved  in  the  model.  Model  geometry  is 
the  same  as  the  one  in  our  previous  study  and  above  figure  is  the  same  as  Fig.  lb  of 
Ref.  [36]. 


the  electrolyte.  As  a  result  the  cell  is  short  circuited  and  drop  in 
open  circuit  voltage  is  experienced.  The  effect  of  current  leaks 
becomes  less  significant  at  higher  current  densities,  because  they 
are  inversely  proportional  to  the  ionic  current  density  [39].  Here, 
we  model  the  electronic  current  leaks  as  boundary  conditions  to 
the  electronic  charge  equation  at  the  electrolyte  interfaces  of  anode 
and  cathode,  which  will  be  explained  later  while  describing  the 
boundary  conditions. 

The  mathematical  model  is  summarized  below.  However,  for  the 
details  of  the  equations  and  the  modeling  framework  that  includes 
a  realistic  model  of  the  actual  fuel  cell  test  environment,  we  refer 
to  our  previous  studies  [36,37]. 


2  A.  Mass  and  momentum  conservation 

Continuity  equation  is  solved  coupled  with  non-isothermal 
momentum  equation  to  account  for  mass  and  momentum  balances 
respectively. 

(„„),£>  (,) 


d_ 

dt 


(^)  +  pu-Vu  =  V- 


-pI+^(Vu  +  (Vu)T) 


^(V-u)I 


+  Pg-(0U  +  F 


(2) 


In  Eq.  (2),  right  hand  side  describes  the  surface  (pressure  and 
shear)  forces,  body  forces,  Darcy’s  term  accounting  for  the  addi¬ 
tional  pressure  drop  in  porous  media  and  any  other  forces  acting 
on  the  fluid,  F.  The  above  equations  are  written  in  a  general  form 
which  can  be  applied  to  both  gas  channels  and  porous  electrodes. 
In  gas  channels,  porosity,  £  is  taken  to  be  unity  while,  permeability, 
I<  is  infinity.  Source  terms  appearing  in  Eqs.  (1)  and  (2)  as  well  as 
the  ones  for  the  other  conservation  equations  are  given  in  Table  1. 


Density  is  calculated  for  a  multicomponent  mixture  by: 

p=i w'Dm  (3) 

j 

where  Xj  is  the  mole  fraction  of  the  species  determined  from  the 
species  balance  and  Mj  is  the  molecular  weight  of  the  species. 


2.2.  Species  conservation 

In  this  work,  Maxwell-Stefan  equations  are  utilized  for  species 
conservation: 


djepwj) 

dt 


+  V- 


n 

pWjU  -  pWj^Djk  (wxk  +  (; Xk  -  wk)^) 
k= 1 


=  Ri 


(4) 


Eq.  (4)  represents  a  set  of  partial  differential  equations  solved 
for  n  -  1  species  in  a  multicomponent  gas  mixture  of  n  species.  The 
mass  fraction  of  the  nth  species  is  calculated  via: 

n 

Ewj = 1  (5) 

j=  i 


In  the  fuel  channel  Eq.  (4)  is  solved  for  hydrogen  and  water, 
whereas  in  the  air  chamber  it  is  solved  for  oxygen  and  water.  Nitro¬ 
gen  is  the  background  (nth)  species  for  both  mixtures.  In  the  form 
of  Maxwell-Stefan  equations  developed  by  Curtiss  and  Bird  [40], 
second  term  inside  the  brackets  on  the  right  hand  side  is  defined 
as  the  diffusion  driving  force.  Derivation  of  the  above  term,  valid 
for  ideal  gases  can  be  found  in  Ref.  [40].  Djk  in  Eq.  (4)  are  the 
Maxwell-Stefan  diffusivities  and  calculated  from  binary  diffusion 
coefficients  for  each  species  pair.  Binary  diffusion  coefficients  are 
calculated  using  Fuller’s  method  [41,42]  which  is  reported  to  match 
the  experimental  data  better  at  high  temperatures  [43]. 

„  0.00143T1-75 

Djk  =  - ?  (6) 

pMf2(vf3 +Vf3) 

in  which  p  is  the  gas  pressure  and  T  is  the  temperature,  in  bars 
and  I<,  respectively.  Djk  is  calculated  in  m2  s-1  from  this  expression. 
Mjk  is  defined  as  the  combined  molecular  weight  of  the  gas  pair 
and  can  be  calculated  as  MJ-fe  =  2[l/MJ- +  1/Mk]_1.  Vj  are  the  special 
diffusion  volumes  of  Fuller  et  al.  [41  ],  which  are  listed  for  SOFC  gases 
in  Ref.  [44].  Binary  diffusivities  are  then  corrected  to  account  for 
mass  transfer  resistance  in  porous  media  via  the  following  relation: 


Dff  =  -D 


jk 


’jk 


(7) 


where  s  is  the  porosity  and  r  is  the  tortuosity  and  listed  in  Table  2. 
Porosity  values  for  unreduced  anode  and  cathode  are  reported  in 
[10]  whereas  tortuosity  values  are  taken  from  [18]. 

In  porous  media,  Eq.  (4)  is  modified  to  implement  the  dusty  gas 
model  (DGM)  due  to  the  impact  of  Knudsen  diffusion  in  small  pores. 
DGM  can  be  formulated  either  using  Fields  Law  or  Maxwell-Stefan 
approach  to  model  the  species  balance  [45].  Using  the  Maxwell- 
Stefan  approach  the  DGM  formulation  can  be  written  as  [46] 


XvNj  -  XjN;  Nj  1  ^ 

J  - LJ.  j - L_  = - Vp, 

rfff  rfff  RT  Pl 

ij  i,Kn 

jj=i 


(8) 


Table  1 

Source  terms. 
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Mass  and  Species 
Momentum 

Energy 

Anode  transfer  current 

Cathode  transfer  current 


Ri=±WM> 

f=uERj 

Q  =  (£max  -  Vcell)ia 

'“=A“exp  (~w) 

ic=Acexp(-E^ 


xh2P 
Xi/2Pref 

X02P 
x02Pref  J 


xh20  P 
xh20  Pref 


sinh 


(  OCaF_ 
^  RT 


Pa 


)  sinh(Sr"c) 


(18) 

(19) 

(20) 
(21) 

(22) 


where  is  the  effective  Knudsen  diffusion  coefficient  which  is 

i,  An 

calculated  from  Mason  and  Malinauskas  [46]  as: 


4  s  1 8RT 

Di,Kn  “  3  ^  y  ^MirP 


(9) 


Knudsen  diffusion  coefficient  is  calculated  for  each  species  sep¬ 
arately.  In  Eq.  (9)  rp  is  the  average  pore  size. 

Eq.  (8)  can  be  cast  into  a  matrix  form  [47] 

[B](n)  =  “rT(Vp)  (10) 

which  is  converted  for  species  diffusion  flux  by  inversion  of  matrix 

[B] 

(N)  =  _Rr[B]~1(Vp)  (11) 

N  is  then  used  in  Eq.  (4)  which  replaces  the  second  term  inside 
the  brackets  on  the  left  hand  side. 


reported  that  neglecting  variations  of  the  gas  stream  properties 
along  the  gas  channels  has  a  significant  impact  on  both  steady  state 
and  transient  results,  therefore  fluid  properties  for  a  single  species 
are  defined  as  functions  of  temperature  using  the  thermodynamic 
property  tables  found  in  [48]. 

In  the  porous  regions  heat  equation  is  formulated  by  volume 
averaging  as  described  in  [49].  Assuming  thermal  equilibrium 
between  the  gas  phase  and  the  solid  phase,  the  volume  averaging 
results  in  the  same  form  of  Eq.  (12),  with  modifications  to  the  tran¬ 
sient  and  conduction  term.  In  the  porous  domains  effective  pCp  and 
thermal  conductivity  are  introduced  in  the  transient  and  conduc¬ 
tion  terms  respectively  such  that  {pCp)eff=s(pCp)f+{\  -  s){pCp)s 
and  keff=  skf+(  1  -  s)ks.  Since  the  thermal  capacity  (pCp)  of  the  solid 
phase  is  much  larger  than  that  of  the  fluid  phase,  contribution  of 
the  fluid  is  neglected  in  the  transient  term. 

The  radiation  heat  transfer  between  the  outer  boundaries  of  the 
cell  (as  well  as  the  alumina  tube)  and  the  furnace  walls  are  modeled 
as  a  boundary  conditions: 


2.3.  Energy  conservation 


n  ■  (qi  -  q2)  =  aabsesb(T4  -  T4all);  ^  =  -/<;VT,  +  p,Cp, ,7,11;  (14) 


The  conservation  of  energy  for  the  entire  domain  is  governed 
by: 

d(p(^r)  +  v  •  (-kVT  +  pCp  7u  +  £h,NDJ)  =  Q  (12) 

j 

where  the  second  term  and  third  term  inside  the  divergence 
account  for  heat  transfer  due  to  bulk  convection  and  diffusion. 
These  terms  disappear  in  solid  regions.  Specific  heat  for  the  fluid 
mixture  is  determined  via  the  expression 

Cp,mix  —  ^  >JCPJ  (13) 

j 

Although  there  are  recommended  expressions  to  estimate  mix¬ 
ture  thermal  conductivity  and  viscosity  [43],  the  values  for  N2  is 
used  for  other  mixture  properties  since  both  the  air  and  fuel  streams 
are  largely  comprised  of  N2.  On  the  other  hand  Ioara  et  al.  [27] 


2.4.  Charge  conservation 


Since  the  characteristic  time  scale  for  charge  transfer  are  much 
smaller  than  those  for  heat  and  mass  transfer  [30],  ionic  and  elec¬ 
tronic  charge  conservation  equations  are  solved  in  steady  state. 


V  •  {oeV<t>e)  =  j 

(15) 

V  •  (OjV$f)  =  -j 

(16) 

Source  terms  in  these  equations  are  the  charge  transfer  rates 
calculated  via  Butler-Volmer  equations  and  are  listed  in  Table  1. 

To  account  for  the  current  leaks  through  the  electrolyte,  the  fol¬ 
lowing  electronic  current  is  given  as  a  flux  boundary  condition  to 
the  electronic  charge  equation  (Eq.  (15))  at  the  anode  electrolyte 
interface: 


n  •  (cTjV^)  _  j/eak,  jleak  ~ 


0eq  —  0e,c 


®i,elc 


Po 

P02,c 


1/4 


Table  2 

Geometrical  parameters  ([10]  unless  stated  otherwise). 


Property 

Value 

Inner  anode  diameter 

0.8  mm 

Anode  thickness 

0.4  mm 

Electrolyte  thickness 

30|jim 

Cathode  thickness 

70|jim 

Anode  length 

10  mm 

Electrolyte  length 

10  mm 

Cathode  length 

7  mm 

Anode  porosity  (after  reduction  of  NiO),  ea 

0.45 

Cathode  porosity,  sc 

0.55 

Anode  tortuosity,  za  [18] 

3.5 

Cathode  tortuosity,  rc  [18] 

3.5 

exp((F//?r)0e,c)  -  1  fl7x 

1-exp  (-(F/RTUeq-(Pe,c))  ^  J 

At  the  cathode  electrolyte  interface,  same  flux  is  given  as  a 
flux  boundary  condition  to  the  electronic  charge  equation.  Eq.  (17) 
derived  by  Riess  [50]  is  the  analytical  solution  to  the  electronic 
charge  transfer  equation  in  the  electrolyte. 

3.  Numerical  method 

Model  equations  are  implemented  in  a  commercial  multiphysics 
software,  COMSOL,  which  uses  finite  element  method  to  discretize 
the  partial  differential  equations  (PDEs).  For  transient  problems 
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Fig.  2.  Current  density  response  of  the  SOFC  to  (i)  single  step  change  in  voltage  from 
0.7  to  0.4  V  at  t=  1  s  (dashed  line),  (ii)  multiple  step  change  in  voltage  from  0.7  to 
0.6  V  at  t=l  s  then  0.6-0.5V  at  t=6s  and  0.5-0.4V  at  t=  11  s  (solid  line).  The  insert 
shows  the  very  initial  time  intervals  of  the  current  response  for  the  single  step  case 
to  show  the  overshoot  phenomena. 

COMSOL  uses  (vertical)  method  of  lines  in  which  it  first  discretizes 
the  PDEs  in  the  space  domain  [51  ].  After  reducing  the  equation  sys¬ 
tem  to  a  set  of  non-linear  ordinary  differential  equations  (ODEs), 
COMSOL  uses  DASPK  (differential-algebraic  solver  preconditioned 
Krylov)  algorithm  to  solve  for  the  new  set  of  equations.  DASPK  algo¬ 
rithm  developed  by  Brown  et  al.  [52]  is  an  implicit  time-stepping 
scheme  which  implies  that  equation  system  of  ODEs  is  solved  in 
each  time  step.  Non-linear  set  of  equations  is  first  solved  with  New¬ 
ton  iteration  and  then  one  of  the  linear  solvers  in  COMSOL’ s  library 
is  employed  in  the  rest  of  the  solution. 

Spatial  discretization  is  carried  out  with  3799  triangular  and 
1940  quadrilateral  mesh  elements.  Use  of  different  element  types  in 
the  same  model  allows  us  to  adjust  the  mesh  to  cut  down  computa¬ 
tional  time  while  ensuring  the  accuracy.  Time  intervals  are  chosen 
as  small  as  10-4  in  the  vicinity  of  the  step  changes  whereas  time 
steps  are  chosen  in  the  order  of  seconds  near  steady  state.  The  solu¬ 
tion  takes  from  5  min  to  1  h  depending  on  the  time  stepping  chosen 
for  specific  analysis. 

4.  Results  and  discussion 

Using  the  model  described  in  the  previous  section,  we  have  sim¬ 
ulated  multiple  transient  scenarios,  first  of  which  is  the  response  to 
a  change  in  operating  voltage.  For  the  following  cases,  the  cell  tem¬ 
perature  is  maintained  at  550  °C  and  the  cell  runs  on  3%  relative 
humidity  fuel  composed  of  5  ccm  hydrogen  and  20  ccm  nitrogen 
while  it  uses  ambient  air  (inside  the  furnace)  on  the  cathode  side. 

4.1  Response  to  a  change  in  cell  voltage 

Fig.  2  shows  the  dynamic  responses  of  the  current  to  a  change 
in  voltage  from  0.7  to  0.4  V  for  two  cases:  (z)  a  single  step  change  in 
voltage  from  0.7  to  0.4  V  at  t=  1  s  (dashed  line);  and  (ii)  a  sequen¬ 
tial  step  change  in  voltage  first  from  0.7  to  0.6  V  at  t  =  1  s,  then  from 
0.6  to  0.5  V  at  t  =  6s  and  then  from  0.5  to  0.4  V  at  t  =  11  s.  In  the 
sequential  case,  second  and  third  steps  are  given  at  times  when  the 
current  density  reaches  98%  of  the  steady-state  values  as  a  result  of 
the  previous  inputs.  For  both  cases  dynamic  response  of  the  cell 
exhibits  a  sudden  jump  in  current  density  at  the  instant  when 
the  step  input  is  given  due  to  the  instantaneous  electrochemical 
reaction  taking  place  in  the  electrodes.  It  should  be  noted  that  the 
double-layer  capacitance  is  not  accounted  in  this  model  as  the  asso¬ 
ciated  time  scales  are  much  smaller  (~10  p,s)  than  those  of  other 


Fig.  3.  Response  of  the  average  cell  temperature  to  the  voltage  changes  described 
in  Fig.  2. 

physical  phenomena  [33].  As  a  result  of  this  momentary  increase  in 
current  density,  reactant  concentration  at  the  reaction  sites  drops 
instantaneously.  Due  to  the  larger  mass  transfer  time  scale  from 
the  channel  to  the  electrodes,  the  reactants  at  the  catalyst  sites  are 
not  replaced  instantly.  The  time  scale  for  mass  transfer  is  propor¬ 
tional  to  the  diffusion  time  (e.g.  82/D ),  is  calculated  to  be  in  the 
order  of  3  x  10-3  s  for  anode  2  x  10-4  s  for  cathode.  Consequently, 
concentration  overpotential  increases  and  hence  the  current  den¬ 
sity  decreases  for  a  period  of  time  which  results  in  an  undershoot 
in  the  dynamic  response  of  the  cell.  This  can  be  seen  in  the  insert  of 
Fig.  2  which  magnifies  the  very  initial  time  intervals  after  the  volt¬ 
age  change  for  single  step  case.  It  is  observed  that  the  undershoot 
lasts  for  0.01  s  which  is  similar  to  the  time  scale  calculated  for  the 
mass  transfer  from  the  gas  channel  to  the  electrodes. 

Following  the  undershoot,  a  transient  regime  of  about  20  s  is 
observed  until  the  system  reaches  a  steady  state.  This  period  is 
mainly  governed  by  the  temperature  and  is  related  to  the  heat  trans¬ 
fer  in  the  fuel  cell.  With  a  sudden  increase  in  current  density,  the 
ohmic  heat  suddenly  increases.  As  a  result,  temperature  distribu¬ 
tion  changes  inside  the  cell,  however  this  process  is  considerably 
slower.  Time  scale  for  heat  transfer  is  in  the  same  order  of  mag¬ 
nitude  with  the  thermal  diffusion  time,  82/a ,  which  is  calculated 
as  40  s.  When  the  single  step  case  is  compared  to  multiple  step 
case,  it  is  seen  that  the  current  response  of  the  latter  is  somewhat 
slower,  although  they  reach  the  same  steady  state.  To  understand 
this  difference,  we  compare  the  temperature  responses.  Increase 
in  average  cell  temperature  with  respect  to  time  is  seen  in  Fig.  3 
for  both  cases.  In  the  single  step  case,  average  temperature  in  the 
positive/electrolyte/negative  (PEN)  rises  from  558  to  567  °C  in  about 
20  s.  In  the  multiple  step  case,  temperature  first  rises  to  560  °C  after 
the  voltage  changes  from  0.7  to  0.6  V  then  it  reaches  to  563  °C  after 
the  second  voltage  change  from  0.6  to  0.5  V  and  then  it  reaches 
567  °C  after  the  voltage  drops  to  0.4  V  at  t=  11  s.  In  this  figure  dotted 
lines  represent  the  temperature  responses  for  the  cell  when  each 
step  change  is  given  as  a  standalone  input  to  the  system.  If  each 
standalone  input  is  compared,  it  is  seen  that  time  scale  for  each 
case  is  identical.  Indeed  for  a  thermal  system,  represented  with 
a  first  order  differential  equation  (Eq.  (12)),  the  time  constant  is 
independent  of  the  step  input  but  is  inversely  related  to  the  effec¬ 
tive  thermal  diffusivity  ( l<lpCp )  of  the  system.  The  magnitude  of  the 
step  only  determines  the  final  steady-state  value.  Hence  time  scale 
for  the  case  when  voltage  drops  from  0.7  to  0.6  V  is  same  as  the 
one  for  single  input  case  when  the  voltage  drops  from  0.7  to  0.4  V 
which  is  about  20  s.  So  at  t=  11  s  when  voltage  is  changed  from  0.5 
to  0.4  V,  it  is  expected  that  the  cell  takes  another  20  s  to  reach  the 
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Fig.  4.  Comparison  of  the  overall  time  scale  of  the  current  density  response  (solid 
line)  with  that  of  temperature  response  (dashed  line)  to  the  single  step  change  in 
cell  voltage  described  in  Fig.  2. 


Fig.  5.  Response  of  the  average  hydrogen  mole  fraction  at  the  anode  channel  outlet 
to  the  voltage  changes  described  in  Fig.  2. 


new  steady  state.  As  a  result,  the  overall  dynamics  of  the  cell  for  the 
sequential  input  case  is  slower  than  that  for  single  input  case. 

Another  important  observation  in  Fig.  3  is  the  rates  of  increase 
in  temperature  for  both  cases.  For  single  input  case  the  increase  in 
temperature  is  sharper  than  that  for  the  sequential  input  case,  since 
for  the  latter  the  heat  generation  in  the  cell  at  each  step  change  is 
much  smaller  than  that  for  the  single  input  case.  Sudden  changes 
in  temperature  have  the  risk  of  thermal  shocks,  which  brittle  SOFC 
materials  may  not  survive.  So  from  the  thermal  stress  perspective, 
multi-step  load  change  may  be  beneficial  although  it  will  delay  the 
cell  performance  steady  state. 

As  mentioned  before  the  overall  dynamics  of  the  cell  is  governed 
by  the  thermal  transients.  Fig.  4  shows  the  correlation  between 
the  transients  in  current  density  and  temperature,  where  the  cur¬ 
rent  density  and  temperature  responses  of  the  cell  to  the  voltage 
change  from  0.7  to  0.4  V  at  t=  1  s  is  plotted.  For  better  comparison, 
current  density  (dashed  line)  and  average  cell  temperature  (solid 
line)  are  normalized  with  the  final  steady-state  values.  Because  the 
dynamics  of  the  current  density  until  1.01  s  is  associated  with  the 
instantaneous  electrochemical  reactions  and  the  transients  of  the 
mass  transfer,  initial  time  intervals  are  not  included  in  the  plot.  As 
seen,  trends  of  both  curves  are  identical  reaching  steady  state  at 
about  20  s.  Note  that  current  density  increase  is  not  linearly  related 
to  the  temperature  increase,  since  parameters  determining  the  cell 
performance  (e.g.  ionic  conductivity  and  exchange  current  density) 
are  exponential  functions  of  temperature. 

Fig.  5  shows  the  average  hydrogen  mole  fraction  at  the  anode 
outlet  as  a  function  of  time  for  single  step  and  multiple  step  changes. 
At  the  instant  of  the  step  change  in  the  voltage,  instantaneous  drop 
in  hydrogen  mole  fraction  is  observed  as  result  of  sudden  jump  in 
the  current  density.  After  this  initial  drop,  mole  fraction  decreases 
gradually  until  the  overall  system  response  becomes  steady  due  to 
the  increasing  current  density. 

Current  density  profiles  at  the  anode-electrolyte  interface  are 
shown  in  Fig.  6  for  the  step  change  in  voltage  from  0.7  to  0.55  V  at 
t=  1  s.  In  Fig.  6a  evolution  of  the  ionic  current  density  profiles  dur¬ 
ing  the  overshoot  is  shown.  Immediately  after  the  change  in  voltage, 
current  density  increases  up  to  1.18  A  cm-2.  After  this  sudden  jump, 
the  profiles  evolve  to  lower  values.  Note  that  the  cathode  coating 
is  shorter  than  the  anode  tube  and  reaction  does  not  take  place  in 
parts  of  the  anode  where  there  is  no  corresponding  cathode  due  to 
limited  finite  ionic  conductivity.  After  the  overshoot,  current  den¬ 
sity  increases  gradually  as  seen  in  Fig.  6b.  After  t=  7  s,  the  change  in 
the  profile  is  not  significant. 


Corresponding  temperature  profiles  at  the  same  interface  are 
seen  in  Fig.  7.  The  effect  of  sudden  jump  in  the  current  density  is  not 
seen  in  the  temperature  response  and  the  profile  at  t=  1.2  s  is  almost 
identical  to  the  initial  steady  state.  With  the  change  in  voltage  from 
0.7  to  0.55  V,  maximum  temperature  in  the  cell  increases  from  561 
to  573  °C.  However,  the  increase  is  more  in  the  first  time  steps  as  the 
maximum  temperature  reaches  565  °C  1.8  s  after  the  voltage  change 


Fig.  6.  Evolution  of  current  density  profiles  along  the  anode-electrolyte  interface 
for  a  change  in  the  cell  voltage  from  0.7  to  0.55  V  at  t=  1  s  a)  during  the  overshoot, 
b)  after  the  overshoot.  Insert  shows  where  the  profiles  are  drawn  in  the  cell. 
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Fig.  7.  Evolution  of  axial  temperature  profiles  along  the  anode-electrolyte  interface 
for  the  change  in  the  cell  voltage  from  0.7  to  0.55  V  at  t=  1  s. 


Fig.  9.  Evolution  of  the  velocity  profile  at  the  inlet  of  the  fuel  channel  during  the 
failure. 


and  then  570  °C  6  s  after  the  voltage  change.  Maximum  temperature 
in  the  cell  is  observed  towards  the  outlet  of  the  mid-point,  since 
convective  flow  in  the  fuel  channel  distorts  the  symmetry  of  the 
temperature  profile.  As  a  result  current  density  profiles  are  also 
distorted. 

For  the  same  input  hydrogen  mole  fraction  profiles  are  seen 
in  Fig.  8  along  the  anode-electrolyte  interface.  Profiles  1-8  show 
the  evolution  of  the  hydrogen  mole  fraction  at  the  interface  until 
1.0051  s.  In  this  period  profiles  change  significantly.  This  change  is 
related  to  the  dynamics  of  the  mass  transfer  from  the  fuel  channel 
to  the  electrode.  Profiles  9,  10  and  11  correspond  to  2,  5  and  40  s 
respectively.  There  is  a  little  difference  between  these  three  pro¬ 
files.  This  slow  change  is  attributed  to  the  slower  dynamics  of  the 
heat  transfer  in  the  cell  which  results  in  a  gradual  increase  in  cur¬ 
rent  generation  hence  a  gradual  decrease  in  hydrogen  mole  fraction. 
Recall  that  regions  from  0  to  1.5  mm  and  8.5  to  10  mm  correspond 
to  inactive  anode  regions. 

4.2.  Failure  in  the  fuel  supply  system 

We  have  also  simulated  a  case  involving  a  malfunction  in  the 
fuel  supply  system  such  as  a  pressure  loss,  due  to  a  failure  in  the 
flow  controllers  or  damage  in  the  fuel  line.  In  this  failure  scenario, 
we  assume  that  fuel  flow  is  shut-down  in  a  very  short  period  of 
time  and  we  investigate  possible  outcomes  of  this  situation. 


Fig.  8.  Evolution  of  hydrogen  mole  fraction  profiles  along  the  anode-electrolyte 
interface  for  the  change  in  the  cell  voltage  from  0.7  to  0.55  V  at  t=  1  s. 


It  is  assumed  that  due  to  the  failure  in  the  fuel  supply  system 
at  t=  1  s,  the  flow  is  completely  stopped  in  0.04s.  Fig.  9  shows 
the  evolution  of  the  velocity  profiles  at  the  fuel  channel  inlet.  As 
inlet  velocity  decreases  less  hydrogen  is  carried  to  the  anode.  As  a 
result  the  current  density  of  the  cell  is  presumed  to  be  decreasing 
abruptly.  Current  density  response  of  a  micro-tubular  SOFC  operat¬ 
ing  at  0.7  V  and  550  °C  to  this  change  in  fuel  supply  at  t  =  1  s  is  seen  in 
Fig.  10.  Recall  that  ceria  based  electrolytes  possess  mixed  conductiv¬ 
ity,  resulting  in  internal  current  leaks.  Consequently,  cell  becomes 
short  circuited  and  the  operating  current  density  is  smaller  than 
the  ionic  current  density.  In  Fig.  10  both  ionic  current  density  (solid 
line)  and  operating  current  density  (dashed  line)  are  plotted.  Even 
though  there  is  still  hydrogen  in  the  fuel  chamber  to  be  used  by 
the  fuel  cell  till  t=  5  s,  the  ionic  current  generated  by  the  cell  is  not 
enough  to  compensate  for  the  current  leaks.  Therefore  in  less  than 
1  s  after  the  failure  in  the  fuel  supply  system  the  fuel  cell  becomes 
incapable  of  generating  useful  (i.e.  external)  current. 

Fig.  11  shows  the  evolution  of  the  hydrogen  mole  fraction  pro¬ 
files  at  the  anode  electrolyte  interface.  Fig.  11a  shows  the  profiles 
till  the  flow  is  completely  stopped.  Fig.  lib  shows  the  profiles  after 
the  total  fuel  flow  rate  becomes  zero.  It  should  be  noted  that  hydro¬ 
gen  and  nitrogen  flow  rates  decrease  with  the  same  rate  thus  the 
hydrogen  mole  fraction  at  the  inlet  stays  constant  until  the  flow 
stops  completely.  While  there  is  flow  at  the  inlet,  hydrogen  is  car- 


Fig.  10.  Ionic  (solid  line)  and  electronic  (dashed  line)  current  density  response  of 
the  SOFC  to  the  cutoff  of  the  fuel  flow  at  t  =  1  s.  In  the  insert  ionic  current  density 
response  is  plotted  in  logarithmic  y-axis. 
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Fig.  11.  Evolution  of  the  hydrogen  mole  fraction  profiles  at  the  anode  electrolyte 
interface  a)  during  the  flow  rate  is  decreasing  to  zero  b)  after  the  fuel  flow  completely 
stops. 


ried  via  both  diffusion  and  convection.  For  25  cm3  min-1  of  total 
fuel  flow  rate  (at  t  =  1  s)  the  maximum  velocity  in  the  center  of  the 
tube  is  3.25  ms-1.  The  effect  of  convective  dominated  process  can 
be  observed  in  hydrogen  mole  fraction  profiles  as  they  are  flatter  in 
the  inactive  anode  regions  and  linear  in  the  active  anode  regions. 
At  t=  1.02  s  mass  transfer  is  still  dominated  by  convection  as  the 
maximum  inlet  velocity  is  still  as  high  as  1.6  m  s-1 .  The  trend  of  the 
hydrogen  mole  fraction  profile  is  very  similar  to  the  one  before  the 
failure.  At  t=  1.03  s  flow  rate  is  not  zero  yet,  however  the  effect  of 
convection  decreases.  This  can  be  seen  in  the  mole  fraction  pro¬ 
file  as  the  curve  becomes  linear  in  the  inactive  anode  regions  and 
parabolic  in  the  active  anode  region. 

After  the  flow  stops  completely  at  t=  1.04  s,  remaining  hydro¬ 
gen  in  the  fuel  channel  is  carried  to  the  anode  purely  by  diffusion. 
Without  the  effects  of  convection,  hydrogen  mole  fraction  profiles 
evolve  to  be  symmetric  as  seen  in  Fig.  lib  and  hydrogen  diffuses 
to  the  anode  from  both  inlet  and  outlet  of  the  tube.  In  only  0.5  s 
hydrogen  mole  fraction  decreases  as  low  as  0.02  at  the  anode  inlet 
but  it  takes  4  s  for  the  cell  to  deplete  nearly  all  the  hydrogen  in  the 
surrounding  gas  channel  completely. 

5.  Conclusion 

A  previously  developed  computational  fluid  dynamics  (CFD) 
based  SOFC  model  [36,37]  is  extended  to  analyze  the  dynamic 
behavior  of  a  micro-tubular  SOFC.  Response  of  the  cell  to  the 
changes  in  voltage  is  investigated.  An  overshoot  is  observed  in  the 


current  density  response  as  a  result  of  the  combined  effect  of  fast 
electrochemical  reaction  and  slower  dynamics  of  the  mass  trans¬ 
fer.  It  is  predicted  that  time-scales  of  a  micro-tubular  SOFC  is  in  the 
order  of  20  s  governed  by  the  dynamics  of  heat  transfer.  Response 
of  the  cell  is  also  compared  for  two  different  voltage  changes:  (i) 
cell  voltage  is  changed  from  0.7  to  0.4  V  as  a  single  step  input, 
(if)  cell  voltage  is  changed  from  0.7  to  0.4  V  as  a  combination  of 
three  sequential  step  inputs.  It  is  observed  that  the  response  of 
the  multiple  input  case  is  slower.  We  also  investigated  the  behav¬ 
ior  of  the  cell  in  case  of  a  failure  in  the  fuel  supply  system.  In  the 
failure  scenario  we  studied,  the  fuel  flow  is  shut-down  abruptly. 
Our  results  show  that  in  case  of  a  failure  in  the  fuel  supply  system 
the  cell  reaction  goes  on  for  another  4  s  till  all  the  hydrogen  in  the 
channels  is  depleted  completely.  Flowever,  the  cell  stops  generat¬ 
ing  output  current  1  s  after  the  failure.  This  is  because  after  then  the 
electrochemical  reaction  is  not  utilized  by  the  external  circuit  due 
to  internal  current  leaks.  Similar  behavior  is  also  expected  during 
normal  shut-down  of  a  ceria  based  SOFC. 
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